#ifndef boundary_H
#define boundary_H

#include <iostream>//for debugging purposes
#include "geometry.h"

namespace ABAN
{
    class EulerBackwardPropellant;

    //typical boundary condition class for fractional step algorithm
    class AbstractBoundaryCondition
    {
        protected:
        SeqType face;
        CellIterator cell;

        public:
        AbstractBoundaryCondition(CellIterator aCell, SeqType aFace):face(aFace),cell(aCell){}
        virtual ~AbstractBoundaryCondition(){}
        virtual void applyEffectOnUMatrixes(EulerBackwardPropellant* propellant,CellContainer* aMesh)=0;
        virtual void applyEffectOnPMatrixes(EulerBackwardPropellant* propellant,CellContainer* aMesh)=0;
        virtual void applyEffectOnCorrectedU(EulerBackwardPropellant* propellant,CellContainer* aMesh, bool updateFromUHat)=0;
        virtual void applyEffectOnPressureGrad(EulerBackwardPropellant* propellant,CellContainer* aMesh)=0;
    };

    //todo: adding external function to dispose boundary pointers...
    class BoundaryContainer:public DatumContainer<AbstractBoundaryCondition *>
    {
    };

    //typical neuman/dirichlet boundary condition class
    class BoundaryCondition:public AbstractBoundaryCondition
    {
        protected:
        struct BoundaryConstants
        {
            FloatType kc,ks;
            BoundaryConstants():kc(0),ks(0){}
            BoundaryConstants(FloatType akc, FloatType aks):kc(akc),ks(aks){}
        };

        FloatType coeff;//a coeffecient that shows direction of A and n

        virtual BoundaryConstants getEqPConstants(EulerBackwardPropellant* propellant,CellContainer* aMesh)=0;
        virtual BoundaryConstants getEqUConstants(SeqType dir, EulerBackwardPropellant* propellant,CellContainer* aMesh)=0;

        public:
        BoundaryCondition(CellIterator aCell, SeqType aFace)
            :AbstractBoundaryCondition(aCell,aFace)
            {coeff=(face<noOfAxis?1.0:-1.0);}

        virtual void applyEffectOnUMatrixes(EulerBackwardPropellant* propellant,CellContainer* aMesh);
        virtual void applyEffectOnPMatrixes(EulerBackwardPropellant* propellant,CellContainer* aMesh);
        virtual void applyEffectOnCorrectedU(EulerBackwardPropellant* propellant,CellContainer* aMesh, bool updateFromUHat);
        virtual void applyEffectOnPressureGrad(EulerBackwardPropellant* propellant,CellContainer* aMesh);
    };

    class PressureZeroGradientBoundary:public BoundaryCondition
    {
        protected:
        virtual BoundaryConstants getEqPConstants(EulerBackwardPropellant* propellant,CellContainer* aMesh);

        public:
        PressureZeroGradientBoundary(CellIterator aCell, SeqType aFace):BoundaryCondition(aCell,aFace){}
    };

    class PressureFreeGradientBoundary:public BoundaryCondition
    {
        protected:
        virtual BoundaryConstants getEqPConstants(EulerBackwardPropellant* propellant,CellContainer* aMesh);

        public:
        PressureFreeGradientBoundary(CellIterator aCell, SeqType aFace):BoundaryCondition(aCell,aFace){}
    };

    class ConstVelBoundary:public PressureZeroGradientBoundary
    {
        protected:
        Vector cteVel;

        virtual BoundaryConstants getEqUConstants(SeqType dir, EulerBackwardPropellant* propellant,CellContainer* aMesh);

        public:
        ConstVelBoundary(CellIterator aCell, SeqType aFace, Vector vel)
            :PressureZeroGradientBoundary(aCell,aFace), cteVel(vel)
            {}
    };

    //outlet boundary with zero gradient velocity
    class SimpleOutlet:public PressureZeroGradientBoundary
    {
        protected:
        virtual BoundaryConstants getEqUConstants(SeqType dir, EulerBackwardPropellant* propellant,CellContainer* aMesh);

        public:
        SimpleOutlet(CellIterator aCell, SeqType aFace)
            :PressureZeroGradientBoundary(aCell,aFace)
            {}
    };

    //convective boundary condition with constant c for each direction
    class SimpleConvectiveOutlet:public PressureZeroGradientBoundary
    {
        protected:
        FloatType c;

        virtual BoundaryConstants getEqUConstants(SeqType dir, EulerBackwardPropellant* propellant,CellContainer* aMesh);

        public:
        SimpleConvectiveOutlet(CellIterator aCell, SeqType aFace, FloatType ac)
            :PressureZeroGradientBoundary(aCell,aFace), c(ac)
            {}
    };

}

#include "prop.h"
namespace ABAN
{
    /*******************************************************/
    /*******************BoundaryCondition*******************/
    /*******************************************************/
    void BoundaryCondition::applyEffectOnPMatrixes(EulerBackwardPropellant* propellant,CellContainer* aMesh)
    {
        //must be extended for brokenfaces
        SeqType cellNo=cell->no;
        SeqType absoluteFace = (face<noOfAxis?face:oppositeSide[face]);

        BoundaryConstants bcc=getEqPConstants(propellant,aMesh);
        BoundaryConstants otherBcc=getEqUConstants(absoluteFace,propellant,aMesh);

        FloatType Ub = coeff * (otherBcc.kc*(*aMesh->velsN[absoluteFace])[cellNo] + otherBcc.ks);
        
        propellant->pCoeffMatrix->setVal(cellNo,cellNo,
            propellant->pCoeffMatrix->getVal(cellNo,cellNo) +
            //diffusive part
            (bcc.kc-1.0) / propellant->density / (cell->lengths[absoluteFace]/2.0) * cell->area[absoluteFace]
            );

        (*propellant->pB)[cellNo] = (*propellant->pB)[cellNo] +
           //diffusive part
             (-bcc.ks) / propellant->density / (cell->lengths[absoluteFace]/2.0) * cell->area[absoluteFace] +
             (+Ub / propellant->dt*cell->area[absoluteFace]);
        
/*        if (cellNo==0)
        {
          printf("face : %d      flux : %f\n",face, Ub*cell->area[absoluteFace]);
          getchar();
        }*/
    }
    //----------------------------------------------------------------------------
    void BoundaryCondition::applyEffectOnUMatrixes(EulerBackwardPropellant* propellant,CellContainer* aMesh)
    {
        SeqType cellNo=cell->no;
        SeqType absoluteFace = (face<noOfAxis?face:oppositeSide[face]);

        for(SeqType dir=dirXP; dir<noOfAxis; dir++)
        {
            BoundaryConstants bcc=getEqUConstants(dir,propellant,aMesh);
            BoundaryConstants otherBcc=getEqUConstants(absoluteFace,propellant,aMesh);

            FloatType Ub = coeff * (otherBcc.kc*(*aMesh->velsN[absoluteFace])[cellNo] + otherBcc.ks);

            propellant->velsCoeffMatrix[dir]->setVal(cellNo,cellNo,
                propellant->velsCoeffMatrix[dir]->getVal(cellNo,cellNo) +
                //diffusive part
                (1.0-bcc.kc) * propellant->nu / (cell->lengths[absoluteFace]/2.0) * cell->area[absoluteFace] +
                //convective part
                bcc.kc * Ub * cell->area[absoluteFace]
                );


// std::cout<<""<<(1.0-bcc.kc) * propellant->nu / (cell->lengths[absoluteFace]/2.0) * cell->area[absoluteFace]


            (*propellant->velsB[dir])[cellNo] = (*propellant->velsB[dir])[cellNo] +
               //diffusive part
               bcc.ks * propellant->nu / (cell->lengths[absoluteFace]/2.0) * cell->area[absoluteFace] +
               //convective part
               (-bcc.ks * Ub * cell->area[absoluteFace]);
        }
    }
    //----------------------------------------------------------------------------
    void BoundaryCondition::applyEffectOnPressureGrad(EulerBackwardPropellant* propellant,CellContainer* aMesh)
    {
        SeqType cellNo=cell->no;
        SeqType dir = (face<noOfAxis?face:oppositeSide[face]);

        BoundaryConstants bcc=getEqPConstants(propellant,aMesh);
        FloatType pressureOnFace=bcc.kc*(*aMesh->p)[cellNo]+bcc.ks;
        FloatType areaF=cell->area[dir];

        (*aMesh->pGradients[dir])[cellNo]+=coeff*areaF*pressureOnFace;
    }
    //----------------------------------------------------------------------------
    void BoundaryCondition::applyEffectOnCorrectedU(EulerBackwardPropellant* propellant,CellContainer* aMesh, bool updateFromUHat)
    {
        FloatType updateLevelCoeff=(updateFromUHat?-1.0:1.0);
        SeqType cellNo=cell->no;
        SeqType dir = (face<noOfAxis?face:oppositeSide[face]);

        BoundaryConstants bcc=getEqPConstants(propellant,aMesh);
        FloatType pressureOnFace=bcc.kc*(*aMesh->p)[cellNo]+bcc.ks;
        FloatType areaF=cell->area[dir];

        (*aMesh->vels[dir])[cellNo]-=updateLevelCoeff*coeff*areaF*pressureOnFace*propellant->dt/propellant->density/cell->volume;
    }
    
    /*******************************************************/
    /*************PressureZeroGradientBoundary**************/
    /*******************************************************/
    BoundaryCondition::BoundaryConstants PressureZeroGradientBoundary::getEqPConstants(EulerBackwardPropellant* propellant,CellContainer* aMesh)
    {
        return BoundaryConstants(1.0,0.0);
    }
    /*******************************************************/
    /*************PressureFreeGradientBoundary**************/
    /*******************************************************/
    BoundaryCondition::BoundaryConstants PressureFreeGradientBoundary::getEqPConstants(EulerBackwardPropellant* propellant,CellContainer* aMesh)
    {
        SeqType absoluteFace(face<noOfAxis?face:oppositeSide[face]);
        FloatType dx(cell->lengths[absoluteFace]/2.0); //dx is distance between cell center to face
        
        return BoundaryConstants(1.0,dx*cell->LSEgradP[absoluteFace]);
    }
    /*******************************************************/
    /********************ConstVelBoundary*******************/
    /*******************************************************/
    BoundaryCondition::BoundaryConstants ConstVelBoundary::getEqUConstants(SeqType dir, EulerBackwardPropellant* propellant,CellContainer* aMesh)
    {
        return BoundaryConstants(0.0,cteVel[dir]);
    }
    /*******************************************************/
    /*********************SimpleOutlet**********************/
    /*******************************************************/
    BoundaryCondition::BoundaryConstants SimpleOutlet::getEqUConstants(SeqType dir,EulerBackwardPropellant* propellant,CellContainer* aMesh)
    {
        return BoundaryConstants(1.0,0.0);
    }
    /*******************************************************/
    /****************SimpleConvectiveOutlet*****************/
    /*******************************************************/
    BoundaryCondition::BoundaryConstants SimpleConvectiveOutlet::getEqUConstants(SeqType dir,EulerBackwardPropellant* propellant,CellContainer* aMesh)
    {
        
/*        SeqType absoluteFace(face<noOfAxis?face:oppositeSide[face]);
        FloatType dx(cell->lengths[absoluteFace]/2.0), //dx is distance between cell center to face
                  dt(propellant->dt);
        return BoundaryConstants(dt*c/(dx+dt*c),dx/(dx+dt*c)*((*aMesh->velsN[dir])[cell->no]+cell->LSEgradV[dir][absoluteFace]*dx));*/
//         return BoundaryConstants(0.0,(dt*c/(dx+dt*c))*(*aMesh->vels[dir])[cell->no]+dx/(dx+dt*c)*((*aMesh->vels[dir])[cell->no]+cell->LSEgradV[dir][absoluteFace]*dx));

        SeqType absoluteFace(face<noOfAxis?face:oppositeSide[face]);
        FloatType kc(0.0),ks(0.0);
        FloatType dpdn=cell->getNormalToFaceGradient(aMesh->p,oppositeSide[face]);
        FloatType Uf=(-coeff)*cell->getFaceValue(aMesh->vels[absoluteFace],oppositeSide[face]);
                
        Uf-=dpdn*propellant->dt/propellant->density;
        
        if (absoluteFace==dir)
            ks=-Uf;//1.5-6.0*cell->center[dirZP]*cell->center[dirZP];
        return BoundaryConstants(kc,ks);
    }
    
}

#endif
